This R code includes the analyses and figures of the SMART manuscript based on the validation cohort of 95 AML patient samples. It includes analyses of Junyan Lu, with extensions and modifications by Peter-Martin Bruch.

1 Load libraries and data

Libraries

library(table1)
## 
## Attache Paket: 'table1'
## Die folgenden Objekte sind maskiert von 'package:base':
## 
##     units, units<-
library(pheatmap)
library(RColorBrewer)
library(ggbeeswarm)
## Lade nötiges Paket: ggplot2
## Warning: Paket 'ggplot2' wurde unter R Version 4.2.3 erstellt
library(grImport)
## Warning: Paket 'grImport' wurde unter R Version 4.2.3 erstellt
## Lade nötiges Paket: grid
## Lade nötiges Paket: XML
## Warning: Paket 'XML' wurde unter R Version 4.2.3 erstellt
library(cowplot)
library(geometry)
## Warning: Paket 'geometry' wurde unter R Version 4.2.3 erstellt
library(tidyverse)
## Warning: Paket 'tidyverse' wurde unter R Version 4.2.3 erstellt
## Warning: Paket 'tibble' wurde unter R Version 4.2.3 erstellt
## Warning: Paket 'readr' wurde unter R Version 4.2.3 erstellt
## Warning: Paket 'dplyr' wurde unter R Version 4.2.3 erstellt
## Warning: Paket 'lubridate' wurde unter R Version 4.2.3 erstellt
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter()    masks stats::filter()
## ✖ dplyr::lag()       masks stats::lag()
## ✖ lubridate::stamp() masks cowplot::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Load processed data

load("../../data/Validation_data.Rdata")

filepath<-"./00_raw_data/"
source("./functions.R")
## Lade nötiges Paket: ggpubr
## Warning: Paket 'ggpubr' wurde unter R Version 4.2.3 erstellt
## 
## Attache Paket: 'ggpubr'
## Das folgende Objekt ist maskiert 'package:cowplot':
## 
##     get_legend
## 
## Attache Paket: 'survminer'
## Das folgende Objekt ist maskiert 'package:survival':
## 
##     myeloma
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
fontsize=18

1.0.1 Functions for patient characteristics table

my.render.cont <- function(x) {
    with(stats.apply.rounding(stats.default(x), digits=3), c("",
        "Median [Min, Max]"=sprintf("%s [%s,  %s]", MEDIAN, MIN, MAX)))
}

my.render.cont_IQR <- function(x) {
    with(stats.apply.rounding(stats.default(x), digits=3), c("",
        "Median [Min, Max]"=sprintf("%s [%s,  %s]", MEDIAN, IQR)))
}

my.render.cat <- function(x) {
    c("", sapply(stats.default(x), function(y) with( 
      y,sprintf("%d (%0.0f %%)", FREQ, PCT))))
}

2 Preprocessing

Cytarabin screened twice with the concentration of 20000, average the values

avgTab <- group_by(screenData, name, Concentration1, Concentration2, patID) %>%
    summarise(normVal = mean(normVal))
## `summarise()` has grouped output by 'name', 'Concentration1', 'Concentration2'.
## You can override using the `.groups` argument.
screenData <- distinct(screenData, Drug1, Concentration1, Drug2, Concentration2, name, wellType, patID) %>%
    left_join(avgTab, by = c("patID","name","Concentration1","Concentration2"))

2.1 Calculate AUC for single drugs

Function for calculating AUC using trapezoidal rule

calcAUC <- function(value, conc) {
    #value is a vector of normalized viability values
    #conc is a vector of concentrations.
    
    valueConc <- data.frame(viab=value, conc=conc)
    valueConc <- valueConc[order(valueConc$conc),]
    areaTotal <- 0
    for (i in seq(length(value)-1)) {
        areaEach <- (valueConc$viab[i] + valueConc$viab[i+1])*log(valueConc$conc[i+1]/valueConc$conc[i])*0.5
        areaTotal <- areaTotal + areaEach
    }
    areaTotal/log(valueConc[nrow(valueConc),]$conc/valueConc[1,]$conc)
}

Calculate AUC

viabTab <- screenData %>%
    filter(wellType == "single") %>%
    mutate(Drug = Drug1,
           conc = Concentration1)

aucTab.single <- group_by(viabTab, patID ,Drug) %>%
    summarise(auc = calcAUC(normVal, conc))
## `summarise()` has grouped output by 'patID'. You can override using the
## `.groups` argument.

2.2 Calculate volumn under surface for combinations

getVolume=function(df) {
  #find triangular tesselation of (x,y) grid
  res=delaunayn(as.matrix(df[,-3]), output.options = TRUE)
  #calulates sum of truncated prism volumes
  sum(mapply(function(triPoints,A) A/3*sum(df[triPoints,3]),
             split.data.frame(res$tri,seq_along(res$areas)),
             res$areas))
}

calcVUS <- function(conc1, conc2, normVal) {
    eachTab <- tibble(conc1 = conc1, conc2=conc2, normVal=normVal)
    df <- eachTab %>%
        mutate(conc1=log10(conc1), conc2 = log10(conc2)) %>%
        mutate(conc1 = ifelse(is.infinite(conc1),0,conc1),
               conc2 = ifelse(is.infinite(conc2), 0, conc2)) %>%
        data.frame()
    dfMax <- df
    dfMax$normVal <- 1
    v1 <- getVolume(df)
    v2 <- getVolume(dfMax)
    return(v1/v2)
}
#get combination
viabTab.combi <- screenData %>%
    filter(wellType == "combi")

#calculate volunm under surface
aucTab.combi <- group_by(viabTab.combi, name, patID) %>%
    summarise(auc = calcVUS(Concentration1, Concentration2, normVal)) %>%
    dplyr::rename(Drug = name)
## `summarise()` has grouped output by 'name'. You can override using the
## `.groups` argument.

2.3 Combine auc of single and combi drugs to one table

aucTab <- bind_rows(aucTab.single, aucTab.combi)

Individual concentrations for visualization

concTab <- screenData %>%
    filter(wellType %in% c("single","combi"))

2.4 Annotate in-vivo response and eln risk

aucTab <- left_join(aucTab, patMeta, by = "patID")
concTab <- left_join(concTab, patMeta, by = "patID")

3 Overview of the dataset

3.1 Pie-chart for ELN risk and in-vivo response

library(webr)
## Warning: Paket 'webr' wurde unter R Version 4.2.3 erstellt
plotTab  <- aucTab %>% distinct(patID, ELN22, group) %>%
    group_by(ELN22) %>%
    mutate(n1 = length(patID)) %>%
    mutate(ELNrisk = sprintf("%s\n(n=%s)",ELN22,n1)) %>%
    group_by(group, ELNrisk) %>%
    summarise(n = length(patID)) %>%
    ungroup() %>%
    dplyr::rename(invivo = group)
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
PieDonutCustom(plotTab, aes(ELNrisk, invivo, count=n, pieLabel = "ELN risk"),
         showPieName = FALSE,
         showRatioDonut = TRUE, showRatioPie = FALSE,
         r0 = getOption("PieDonut.r0", 0.2),
         r1 = getOption("PieDonut.r1", 0.7),
         r2 = getOption("PieDonut.r2", 1), piLabel = "ELN risk", donutLabel = "in vivo response",
         donutLabelSize = 4)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

aucTab %>% distinct(patID, ELN22, group) %>% 
  group_by(ELN22, group) %>% 
  count()
## # A tibble: 6 × 3
## # Groups:   ELN22, group [6]
##   ELN22        group             n
##   <chr>        <fct>         <int>
## 1 adverse      Non-Responder    28
## 2 adverse      Responder         8
## 3 favorable    Non-Responder     6
## 4 favorable    Responder        10
## 5 intermediate Non-Responder    14
## 6 intermediate Responder        29
writexl::write_xlsx(rownames_to_column(as.data.frame(plotTab),var="rowname"), path=paste0("../../","docs/Source data/","Fig4A",".xlsx"))

3.2 Patient characteristics table

patMeta_Validation<-left_join(patMeta_Validation, patMeta, by="patID")

label(patMeta_Validation$AGE) <-"Age (years)"
label(patMeta_Validation$diagnosis) <-"Diagnosis"
label(patMeta_Validation$trtstr) <-"Prescribed treatment"
label(patMeta_Validation$SEX) <-"Sex"
label(patMeta_Validation$BMB) <-"Tumor cell infiltration"
label(patMeta_Validation$group) <-"In vivo response group"
label(patMeta_Validation$ELN22) <-"ELN-22 risk group"

# table 1 printing
table1(~ SEX + AGE + diagnosis + BMB + trtstr + group + ELN22,  data=patMeta_Validation, overall="Total",  render.continuous=my.render.cont, render.categorical=my.render.cat)
Total
(N=95)
Sex
female 47 (49 %)
male 48 (51 %)
Age (years)
Median [Min, Max] 59.0 [18.0, 84.0]
Diagnosis
AML 95 (100 %)
Tumor cell infiltration
Median [Min, Max] 76.0 [50.0, 99.0]
Prescribed treatment
Cytarabin + Daunorubicin 95 (100 %)
In vivo response group
Non-Responder 48 (51 %)
Responder 47 (49 %)
ELN-22 risk group
adverse 36 (38 %)
favorable 16 (17 %)
intermediate 43 (45 %)
Supp_Table_Val_Cohort<-patMeta_Validation %>% 
  rename(PatientID=patID, Age=AGE, Sex=SEX, Treatment=trtstr, `Tumor cell infiltration`=BMB, Diagnosis=diagnosis, `In vivo response group`=group, `ELN-22 risk group`=ELN22 ) %>% 
  select(-EFSSTAT, -EFSTM, -EFSTMU) %>% 
  select(PatientID, Age, Sex, Treatment, `Tumor cell infiltration`, Diagnosis, `In vivo response group`, `ELN-22 risk group`) 

writexl::write_xlsx(Supp_Table_Val_Cohort,path="../../docs/table_S3_Validation_patients.xlsx")

3.3 Drug-drug correlation heatmap (Supplementary figure)

viabMat <- aucTab %>% select(patID, Drug, auc) %>%
    pivot_wider(names_from = patID, values_from = auc) %>%
    column_to_rownames("Drug") %>% as.matrix()
corMat <- cor(t(viabMat))
pheatmap::pheatmap(corMat, color = colorRampPalette(c("#003DA5", "white", "#A6093D"))(100), 
         filename=paste0("../../docs/Figure_ED2.pdf"),
                   clustering_method = "ward.D2", main = "Drug-Drug correlation heatmap")

writexl::write_xlsx(rownames_to_column(as.data.frame(corMat),var="rowname"), path=paste0("../../","docs/Source data/","ED2",".xlsx"))

3.4 Drug response heatmap

colAnno <-  aucTab %>% distinct(patID, ELN22, group) %>%
    mutate(ELN22 = factor(ELN22, levels = c("adverse","intermediate","favorable")),
           group = factor(group, levels = c("Non-Responder","Responder"))) %>%
    arrange(ELN22, group) %>%
    dplyr::rename(`ELN risk` = ELN22, `in vivo response`= group) %>%
    column_to_rownames("patID") 

annoColor <- list(`ELN risk` = c(adverse = "#E18727FF", intermediate = "#6F99ADFF", favorable = "#20854EFF"),
                  `in vivo response` = c(`Non-Responder` = "orangered3", Responder = "deepskyblue4"))
viabMat <- viabMat[,rownames(colAnno)]
pheatmap(viabMat, annotation_col = colAnno,
         cluster_cols = FALSE, annotation_colors = annoColor,
         # filename=paste0("../../docs/validation_drugHeatmap.pdf"),
         show_colnames = FALSE,
         treeheight_row = 20)

With hierarchical clustering

viabMat <- viabMat[,rownames(colAnno)]
viabMat.scaled <-  mscale(viabMat, censor = 4)
p <- pheatmap(viabMat.scaled, annotation_col = colAnno,
         cluster_cols = TRUE, annotation_colors = annoColor,
         color=colorRampPalette(c("deepskyblue3", "white", "orangered4"))(100), # colorRampPalette(c(blue2, "white", red2))(100)
         clustering_method = "ward.D2",
         main = "Hierarchical clustering of ex vivo responses in AML validation cohort", 
         show_colnames = FALSE,
         treeheight_row = 20, silent = TRUE)
cowplot::plot_grid(p$gtable)

# ggsave("../../docs/HC_heatmap_Validation.png", plot = p$gtable, height = 8, width =15)
ggsave("../../docs/Figure_ED4.pdf", plot = p$gtable, height = 8, width =15)


writexl::write_xlsx(rownames_to_column(as.data.frame(viabMat.scaled),var="rowname"), path=paste0("../../","docs/Source data/","ED4",".xlsx"))

4 Test correlations between ex-vivo response and in-vivo responses

4.1 Univariate test (without blocking for ELN risk)

testTab <- aucTab
resTab.uni <- group_by(testTab, Drug) %>% nest() %>%
    mutate(m = map(data, ~t.test(auc~group,., var.equal=TRUE))) %>%
    mutate(res = map(m, broom::tidy)) %>%
    unnest(res) %>%
    select(Drug, estimate, p.value) %>%
    arrange(p.value) %>% ungroup() %>%
    mutate(p.adj = p.adjust(p.value, method ="BH"))
sigDrugs <- filter(resTab.uni, p.adj <= 0.1)$Drug

4.1.1 Volcano plot

fdrCut=0.1

plotTab <- resTab.uni %>%
    mutate(dir = ifelse(estimate > 0, "more resistant", "more sensitive")) %>%
    mutate(dir = ifelse(p.adj > fdrCut, "n.s.", dir))
pCut <- -log10(min(filter(resTab.uni, p.adj > fdrCut)$p.value))

pVolcano <- ggplot(plotTab, aes(x=estimate, y=-log10(p.value))) +
    geom_point(aes(fill = dir), shape = 21) + 
    scale_fill_manual(values = c(`more resistant` = "orangered3", `more sensitive` = "deepskyblue3", n.s. = "grey80"), name = "") +
    geom_hline(yintercept = pCut, color = "grey50", linetype = "dashed") + 
    ggrepel::geom_text_repel(data = filter(plotTab, p.adj <= fdrCut),aes(label = gsub("_"," + ",Drug)), size=5) +
    ylab(expression(-log[10]*'('*italic(P)~value*')')) + xlab("Mean viability difference") +
    ggtitle(bquote(italic("In vivo")~"response"~"~"~italic("ex vivo")~"response")) +
    theme_classic()+
    theme(plot.title = element_text(hjust=0.5, face = "bold", size=18),
          axis.title = element_text(size=16),
          axis.text = element_text(size=16),
          legend.text = element_text(size=14),
          legend.background = element_rect(fill = NA),
          axis.line = element_blank(),
          legend.position = c(0.75,0.2),
          plot.margin = margin(0,20,0,20),
          panel.border = element_rect(size=1.5, fill = NA)
          )
        
pVolcano

4.1.2 Boxplot

resTab.sig <- filter(resTab.uni, p.adj <= 0.1)
plotList_box <- lapply(seq(nrow(resTab.sig)) ,function(i){
  rec <- resTab.sig[i,]
  
  #get the df
  df_plot <- filter(aucTab, Drug == rec$Drug)
  pVal <- formatC(rec$p.value, digits = 2)
  drugName <- rec$Drug
  #make plot
  p <- ggplot(df_plot, aes(x=group, y=auc)) + 
    geom_boxplot(aes(fill = group), width = 0.8, size=0.8, alpha = 0.85, 
                 color = "black", outlier.shape = NA) +
    scale_fill_manual(values = c("Non-Responder"="orangered3","Responder"="deepskyblue4")) +
    geom_beeswarm(cex=3, dodge.width=0.8, shape=21, 
                  fill="slategray4", size=4, alpha= 0.7) +
    geom_hline(yintercept = 1, linetype="22", col="gray30", size=0.7) +
    coord_cartesian(ylim = c(0,1.5)) +
    #ggtitle(sprintf("%s (P = %s)",rec$Drug, ) +
    ggtitle(bquote(.(drugName)~"("*italic("P")~"="~.(pVal)~")")) +
    #xlab(expression(paste(italic("in vivo")," response"))) +
    xlab("") +
    ylab("Viability") +
    guides(color="none", fill="none") +
    theme_classic()+
    theme(plot.title = element_text(hjust=0, face = "bold", size=18),
          axis.title = element_text(size=16),
          axis.text = element_text(size=16),
          legend.text = element_text(size=14),
          legend.background = element_rect(fill = NA),
          axis.line = element_blank(),
          legend.position = c(0.75,0.2),
          plot.margin = margin(0,20,0,20),
          panel.border = element_rect(size=1.5, fill = NA)
          )

  return( p )
})
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
names(plotList_box) <- resTab.sig$Drug
pBox <- cowplot::plot_grid(plotlist = plotList_box, ncol=2)
pBox

writexl::write_xlsx(rownames_to_column(as.data.frame(aucTab),var="rowname"), path=paste0("../../","docs/Source data/","ED8A",".xlsx"))

4.2 Dose-response curves for selected drugs

example_drugs <- sigDrugs[sigDrugs!= "Daunorubicin_Cytarabin"]

plotList_DResp <- lapply(example_drugs ,function(drugname){
    
  plotTab <- filter(concTab, Drug1 == drugname, Drug2 == "DMSO")
  
  ciTab <- group_by(plotTab, group, Concentration1) %>%
      summarise(sdVal = sd(normVal), n=length(patID), normVal = mean(normVal)) %>%
      mutate(se = sdVal/sqrt(n), ci = 1.96*se)
  
  #plot values: real with dots, predicted with line
  ggplot(plotTab,aes(x= Concentration1, y=normVal, color=group, fill=group)) + 
    geom_point(shape = 21, size=3.5, alpha = 0.2, color="black") +
    geom_errorbar(data = ciTab, aes(x=Concentration1, y=normVal, ymax = normVal + ci, ymin=normVal-ci), width=0.1, size=1) +
    scale_color_manual(values = c("Responder" = "deepskyblue4", "Non-Responder"="orangered3")) +
    scale_fill_manual(values = c("Responder" = "deepskyblue4", "Non-Responder"="orangered3")) +
    coord_cartesian(ylim = c(0,1.5)) +
    geom_smooth(data = ciTab, method = "fitIC50", se = FALSE) +
    scale_x_continuous(trans='log10') +
    labs(title =  paste0("",drugname) )+
    xlab(expression("Concentration"~(mu*M))) + ylab("Viability") +
    guides(fill="none", color="none") +
    theme_classic()+
    theme(plot.title = element_text(hjust=0, face = "bold", size=18),
          axis.title = element_text(size=16),
          axis.text = element_text(size=16),
          legend.text = element_text(size=14),
          legend.background = element_rect(fill = NA),
          axis.line = element_blank(),
          legend.position = c(0.75,0.2),
          plot.margin = margin(0,20,0,20),
          panel.border = element_rect(size=1.5, fill = NA)
          )
})
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
pDose <- cowplot::plot_grid(plotlist = plotList_DResp, ncol=2)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
pDose

writexl::write_xlsx(rownames_to_column(as.data.frame(concTab),var="rowname"), path=paste0("../../","docs/Source data/","ED8B",".xlsx"))

4.3 Multi-variate test (with blocking for ELN risk)

testTab <- filter(aucTab, Drug %in% sigDrugs)
resTab.multi <- group_by(testTab, Drug) %>% nest() %>%
    mutate(m = map(data, ~car::Anova(lm(auc~group + ELN22,.)))) %>%
    mutate(res = map(m, broom::tidy)) %>%
    unnest(res) %>% filter(term == "group") %>%
    select(Drug, p.value) %>%
    arrange(p.value) %>% ungroup() %>%
    mutate(p.adj = p.adjust(p.value, method ="BH"))
head(resTab.multi)
## # A tibble: 6 × 3
##   Drug                    p.value   p.adj
##   <chr>                     <dbl>   <dbl>
## 1 Vindesine              0.000945 0.00330
## 2 Fludarabine            0.00152  0.00330
## 3 Vincristine            0.00163  0.00330
## 4 Cladribine             0.00189  0.00330
## 5 NSC 652287             0.0133   0.0186 
## 6 Daunorubicin_Cytarabin 0.0465   0.0542

4.4 Test associations between ex-vivo and in-vivo responses, within each ELN group

testTab <- filter(aucTab, Drug %in% resTab.sig$Drug)

resTab <-  testTab %>%
    group_by(Drug, ELN22) %>% nest() %>%
    mutate(m = map(data, ~t.test(auc~group,.,var.equal=TRUE))) %>%
    mutate(res = map(m, broom::tidy)) %>%
    unnest(res) %>% 
    select(Drug, p.value, estimate) %>%
    arrange(p.value) %>% ungroup() %>%
    mutate(p.adj = p.adjust(p.value, method ="BH"))
## Adding missing grouping variables: `ELN22`
head(resTab)
## # A tibble: 6 × 5
##   ELN22   Drug                   p.value estimate  p.adj
##   <chr>   <chr>                    <dbl>    <dbl>  <dbl>
## 1 adverse Vindesine              0.00331    0.249 0.0331
## 2 adverse Fludarabine            0.00417    0.136 0.0331
## 3 adverse Vincristine            0.00473    0.208 0.0331
## 4 adverse Cladribine             0.00663    0.233 0.0348
## 5 adverse Daunorubicin_Cytarabin 0.0138     0.190 0.0579
## 6 adverse NSC 652287             0.0443     0.183 0.155
writexl::write_xlsx(rownames_to_column(as.data.frame(testTab),var="rowname"), path=paste0("../../","docs/Source data/","ED9B",".xlsx"))


writexl::write_xlsx(rownames_to_column(as.data.frame(testTab),var="rowname"), path=paste0("../../","docs/Source data/","Fig4C",".xlsx"))

There are significant p-values, but p.adj is borderline significant

4.4.0.1 Plot p-value heatmap

plotTab <- resTab %>% mutate(pSign = -log10(p.value) * sign(estimate)) %>%
    mutate(ifSig = case_when(p.adj <=0.1 ~ "*",
                             TRUE ~ ""))

pPheatmap <- ggplot(plotTab, aes(x=factor(Drug), y=factor(ELN22, levels = c("adverse", "intermediate", "favorable")), fill = pSign)) +
    geom_tile() + 
    geom_text(aes(label = ifSig), col = "black", size=8, face="bold") +
    scale_fill_gradient2(low = "blue", mid = "white",high = "red", name = bquote("-log10("*italic("P")*")")) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
    xlab("") + ylab("ELN risk group") +
    theme(axis.title = element_text(size=18),
          axis.text = element_text(size=16),
          plot.margin = margin(5,20,5,20))
## Warning in geom_text(aes(label = ifSig), col = "black", size = 8, face =
## "bold"): Ignoring unknown parameters: `face`
pPheatmap 

4.4.0.2 Boxplot of significant associations

resTab.sig <- filter(resTab, p.adj <= 0.1)
pList <- lapply(unique(resTab.sig$Drug), function(dd) {
  
    plotTab <- filter(testTab, Drug == dd) %>%
        mutate(ELN22 = factor(ELN22, levels = rev(c("favorable","intermediate","adverse"))))
    
    pTab <- filter(resTab, Drug == dd) %>%
        mutate(pText = sprintf("italic(P)~'='~%s",formatC(p.value, digits=2)))
    
    ggplot(plotTab, aes(x=group, y=auc))  +
        geom_boxplot(aes(fill = group), width = 0.8, size=0.8, alpha = 0.85, color = "black", outlier.shape = NA) +
        scale_fill_manual(values = c("Non-Responder"="orangered3","Responder"="deepskyblue4")) +
        geom_beeswarm(cex=3, dodge.width=0.8, shape=21, 
                      fill="slategray4", size=4, alpha= 0.7) +
        geom_text(data = pTab, aes(label = pText), x=2.1,y=1.3, vjust=1, parse= TRUE,size=5)+
        geom_hline(yintercept = 1, linetype="22", col="gray30", size=0.7) +
        coord_cartesian(ylim = c(0,1.5)) +
        ggtitle(gsub("_"," + ",dd)) +
        xlab(expression(paste(italic("in vivo")," response"))) +
        ylab("Viability (AUC)") +
        guides(color="none", fill="none") +
        theme_classic()+
        facet_wrap(~factor(ELN22, levels = rev(c("favorable","intermediate","adverse")))) +
        theme(plot.title = element_text(hjust=0, face = "bold", size=18),
              axis.title = element_text(size=16),
              axis.text = element_text(size=16),
              axis.text.x =element_text(angle = 45, hjust = 1),
              legend.text = element_text(size=14),
              legend.background = element_rect(fill = NA),
              axis.line = element_blank(),
              legend.position = c(0.75,0.2),
              plot.margin = margin(0,20,0,20),
              panel.border = element_rect(size=1.5, fill = NA),
              strip.text.x = element_text(size = 16)
              )

})
names(pList) <- unique(resTab.sig$Drug)
#cowplot::plot_grid(plotlist = pList, ncol=1)
mainDrug <- c("Vincristine","Daunorubicin + Cytarabin", "Daunorubicin_Cytarabin")

pList$Daunorubicin_Cytarabin$labels$y<-"Viability (VUC)"
pList$Daunorubicin_Cytarabin$labels$title<-"Daunorubicin + Cytarabine"


pBoxELN_main <- plot_grid(plotlist = pList[names(pList) %in% mainDrug],ncol=1)
pBoxELN_supp <- plot_grid(plotlist= pList[!names(pList) %in% mainDrug],ncol=1)


pBoxELN_main

4.4.1 Dose-response curves for selected drugs per elnGroup

example_drugs <- sigDrugs[sigDrugs!= "Daunorubicin_Cytarabin"]

plotList_DResp <- lapply(example_drugs ,function(drugname){
    
  plotTab <- filter(concTab, Drug1 == drugname, Drug2 == "DMSO") %>%
      mutate(ELN22 = factor(ELN22, levels=c("adverse","intermediate","favorable")))
  
  meanTab <- group_by(plotTab, Drug1, ELN22, group, Concentration1) %>%
      summarise(normVal = mean(normVal, na.rm=TRUE))
  
  #plot values: real with dots, predicted with line
  ggplot(plotTab,aes(x= Concentration1, y=normVal, color=group, fill=group)) + 
    geom_point(shape = 21, size=3.5, alpha = 0.5, color="black") +
    scale_color_manual(values = c("Responder" = "deepskyblue4", "Non-Responder"="orangered4")) +
    scale_fill_manual(values = c("Responder" = "deepskyblue4", "Non-Responder"="orangered4")) +
    coord_cartesian(ylim = c(0,1.5)) +
    geom_smooth(data=meanTab, method = "fitIC50", se = FALSE) +
    scale_x_continuous(trans='log10') +
    labs(title =  paste0("",drugname) )+
    xlab(expression("Concentration"~(mu*M))) + ylab("Viability") +
      facet_wrap(~ELN22) +
    guides(fill="none", color="none") +
    theme_classic()+
    theme(plot.title=element_text(size=fontsize+4, hjust = 0, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
        axis.text.x =element_text(size=fontsize),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_blank(), 
        panel.border = element_rect(size=1.5, fill = NA),
        axis.ticks = element_line(size=1.5),
        legend.text = element_text(size=fontsize),
        legend.title = element_text(size=fontsize+2))
})
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
cowplot::plot_grid(plotlist = plotList_DResp, ncol=1)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'

5 Associate drug responses with clinical outcomes

5.1 Prepare dataset

Outcome data

# clinicTab <- readxl::read_xlsx(paste0("../../data/Validation_cohort_survival.xlsx")) %>%
clinicTab <- patMeta_Validation %>% 
    select(patID, EFSSTAT, EFSTM) %>%
    mutate(across(-patID, as.numeric)) %>%
    dplyr::rename(patID = patID) %>%
    mutate(patID = as.character(patID)) %>%
    distinct(patID, .keep_all = TRUE)

#censor survival data at five years
clinicTab <- mutate(clinicTab, EFSSTAT = ifelse(EFSTM>=60,0,EFSSTAT), EFSTM = ifelse(EFSTM > 60, 60, EFSTM)) %>%
    mutate(EFSTM = EFSTM/12)

ELN data

elnTab <- patMeta %>%
    distinct(patID, ELN22) %>%
    mutate(ELN22 = factor(ELN22, levels = c("favorable","intermediate","adverse")))

survTab <- left_join(elnTab, clinicTab, by = "patID")

Drug response data: use AUC for single drugs and individual concentrations for combinatorial drugs

drugTab <- aucTab %>% 
    filter(Drug %in% sigDrugs)

5.2 Test whether individual drug response associated with survival using univariate model

testTab <- survTab %>% pivot_longer(-c(patID,ELN22)) %>%
    mutate(stat = ifelse(str_detect(name, "STAT"),"endpoint","time"),
           outcome = str_remove(name, "STAT|TM")) %>%
    select(-name) %>% filter(!is.na(value)) %>%
    pivot_wider(names_from = stat, values_from = value) %>%
    left_join(select(drugTab, patID, Drug, auc))
## Joining with `by = join_by(patID)`
resTab <- group_by(testTab, outcome, Drug) %>% nest() %>%
    mutate(m = map(data, ~com(.$auc, .$time,.$endpoint))) %>%
    unnest(m) %>%
    arrange(p) %>%
    select(outcome, Drug, p, HR, lower, higher) %>%
    ungroup() %>%
    mutate(p.adj = p.adjust(p, method ="BH"))

resTab.sig <- filter(resTab, p<0.05)
resTab.sig
## # A tibble: 5 × 7
##   outcome Drug                         p    HR lower higher  p.adj
##   <chr>   <chr>                    <dbl> <dbl> <dbl>  <dbl>  <dbl>
## 1 EFS     Fludarabine            0.00728 13.7   2.03   93.1 0.0257
## 2 EFS     Cladribine             0.0104   4.71  1.44   15.4 0.0257
## 3 EFS     Vincristine            0.0167   6.75  1.41   32.3 0.0257
## 4 EFS     Vindesine              0.0168   5.50  1.36   22.2 0.0257
## 5 EFS     Daunorubicin_Cytarabin 0.0183   5.28  1.32   21.1 0.0257
writexl::write_xlsx(rownames_to_column(as.data.frame(testTab),var="rowname"), path=paste0("../../","docs/Source data/","Fig4D",".xlsx"))

pList <- lapply(seq(nrow(resTab.sig)), function(i) {
    rec <- resTab.sig[i,]
    plotTab <- filter(testTab, outcome == rec$outcome, Drug == rec$Drug) 
    km(plotTab$auc, plotTab$time, plotTab$endpoint,
       gsub("_"," + ",rec$Drug), pval = rec$p, stat = "maxstat", showTable = FALSE)
})
names(pList) <- unique(resTab.sig$Drug)
pKM_main <- pList[names(pList) %in% mainDrug]
pKM_supp <- pList[!names(pList) %in% mainDrug]
pList_withtable <- lapply(seq(nrow(resTab.sig)), function(i) {
    rec <- resTab.sig[i,]
    plotTab <- filter(testTab, outcome == rec$outcome, Drug == rec$Drug)
    km(plotTab$auc, plotTab$time, plotTab$endpoint,
       rec$Drug, pval = rec$p, stat = "maxstat", showTable = TRUE)
})
names(pList_withtable) <- unique(resTab.sig$Drug)
pList_withtable[names(pList_withtable) %in% mainDrug]
## $Vincristine

## 
## $Daunorubicin_Cytarabin

pList_withtable[!names(pList_withtable) %in% mainDrug]
## $Fludarabine

## 
## $Cladribine

## 
## $Vindesine

Supplementary Figure

cowplot::plot_grid(plotlist =pKM_supp, ncol=2)
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet

## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet

## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet

Main figure panel

pLegend <- get_legend(pKM_main$Vincristine)
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
pKM_main$Daunorubicin_Cytarabin$labels$title<-"Daunorubicine + Cytarabine"

pKM_main <- plot_grid(plot_grid(pKM_main$Vincristine + theme(legend.position = "none"),
                      pKM_main$Daunorubicin_Cytarabin + theme(legend.position = "none"), ncol = 2),
                      pLegend, ncol=1 ,rel_heights = c(1,0.1))
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet

## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet

5.3 Test whether individual drug response associated with survival using univariate model per ELN risk group

resTab <- group_by(testTab, outcome, Drug, ELN22) %>% nest() %>%
    mutate(m = map(data, ~com(.$auc, .$time,.$endpoint))) %>%
    unnest(m) %>%
    arrange(p) %>%
    select(outcome, Drug, p, HR, lower, higher) %>%
    ungroup() %>%
    mutate(p.adj = p.adjust(p, method ="BH"))
## Adding missing grouping variables: `ELN22`
resTab.sig <- filter(resTab, p<0.05)
resTab.sig
## # A tibble: 4 × 8
##   ELN22   outcome Drug                         p    HR lower higher  p.adj
##   <fct>   <chr>   <chr>                    <dbl> <dbl> <dbl>  <dbl>  <dbl>
## 1 adverse EFS     Cladribine             0.00355 17.1   2.54  115.  0.0541
## 2 adverse EFS     Fludarabine            0.00515 86.9   3.81 1984.  0.0541
## 3 adverse EFS     Vincristine            0.0229  14.4   1.45  143.  0.160 
## 4 adverse EFS     Daunorubicin_Cytarabin 0.0326   7.32  1.18   45.5 0.171
writexl::write_xlsx(rownames_to_column(as.data.frame(testTab),var="rowname"), path=paste0("../../","docs/Source data/","ED10",".xlsx"))

pList <- lapply(seq(nrow(resTab.sig)), function(i) {
    rec <- resTab.sig[i,]
    leg <- NULL
    pELN <- lapply(c("adverse","intermediate","favorable"), function(eln) {
        eachRec <- filter(resTab, ELN22 == eln, outcome == rec$outcome, Drug == rec$Drug)
        plotTab <- filter(testTab, outcome == rec$outcome, Drug == rec$Drug, ELN22== eln)
        kmPlot <- km(plotTab$auc, plotTab$time, plotTab$endpoint,
           paste0(eln," group"), pval = eachRec$p, stat = "maxstat", showTable = TRUE)
        leg <<- cowplot::get_legend(kmPlot)
        kmPlot <- kmPlot + theme(legend.position = "none")
        kmPlot
        
    })
    
    gTitle <- ggplot() + 
        labs(title = rec$Drug) + 
        theme_void() + theme(plot.title = element_text(hjust = 0.5, face="bold", size= 15))
    cowplot::plot_grid(gTitle,cowplot::plot_grid(plotlist = pELN, nrow=1), leg, ncol=1, rel_heights = c(0.1,1,0.2))
})
cowplot::plot_grid(plotlist = pList, ncol=1)

5.4 Test outcome prediction based on drug response plus ELN risk using multi-variate model

resTab <- group_by(testTab, outcome, Drug) %>% nest() %>%
    mutate(m = map(data, ~coxph(
        Surv(time, endpoint) ~
            auc + ELN22,
        data = .))) %>%
    mutate(res=map(m, broom::tidy)) %>%
    unnest(res) %>%
    filter(term == "auc") %>%
    arrange(p.value) %>%
    select(outcome, Drug, p.value)

resTab.sig <- filter(resTab, p.value<0.05)
resTab.sig
## # A tibble: 5 × 3
## # Groups:   outcome, Drug [5]
##   outcome Drug                   p.value
##   <chr>   <chr>                    <dbl>
## 1 EFS     Fludarabine            0.00537
## 2 EFS     Cladribine             0.0185 
## 3 EFS     Vincristine            0.0211 
## 4 EFS     Daunorubicin_Cytarabin 0.0358 
## 5 EFS     Vindesine              0.0376

5.5 Plot significant associations with forest plots

plotHazard <- function(survRes, title = "") {
    sumTab <- summary(survRes)$coefficients
    confTab <- summary(survRes)$conf.int
    #correct feature name
    nameOri <- rownames(sumTab)
    plotTab <- tibble(feature = c(rownames(sumTab), "reference (favorable)"),
                      HR = c(sumTab[,2],1),
                      p = c(sumTab[,5],1),
                      Upper = c(confTab[,4],1),
                      Lower = c(confTab[,3],1)) %>%
        mutate(Upper = ifelse(Upper > 10, 10, Upper)) %>%
        mutate(feature = case_when(feature == "auc" ~ "ex vivo resistance",
                                   feature == "ELN22intermediate" ~ "intermediate",
                                   feature == "ELN22adverse" ~ "adverse",
                                   TRUE ~ feature)) %>%
        mutate(pLabel = ifelse(!str_detect(feature, "reference"), 
                               sprintf("italic(P)~'='~'%s'", formatNum(p, digits = 1)),
                               "")) %>%
        mutate(feature = factor(feature, levels = c("reference (favorable)","intermediate","adverse","ex vivo resistance")))
    
    ggplot(plotTab, aes(x=feature, y = HR)) +
        geom_hline(yintercept = 1, linetype = "dotted", color = "grey50") +
        geom_point(position = position_dodge(width=0.8), size=3, color = "black") +
        geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.0, size=0.5,color = "grey20") +
        geom_text(position = position_nudge(x = 0.3),
                  aes(y = HR, label =  pLabel),
                  color = "black", size =5, parse = TRUE) +
        expand_limits(y=c(-0.5,0))+
        ggtitle(title) +
        ylab("Hazard ratio") +
        coord_flip() +
        theme_full +
        theme(legend.position = "none", axis.title.y = element_blank())
}
pList <- lapply(seq(nrow(resTab.sig)), function(i) {
    rec <- resTab.sig[i,]
    eachTab <- filter(testTab, outcome == rec$outcome, Drug == rec$Drug) %>%
        select(time, endpoint, auc, ELN22) %>%
        mutate(auc = (auc-mean(auc))/(2*sd(auc)))
    survRes <- runCox(eachTab)
    plotHazard(survRes, paste0(rec$outcome, " ~ ", rec$Drug))
})

names(pList) <- gsub("_"," + ", resTab.sig$Drug)

pList$`Daunorubicin + Cytarabin`$labels$title<-"EFS ~ Daunorubicin + Cytarabine"


writexl::write_xlsx(rownames_to_column(as.data.frame(testTab),var="rowname"), path=paste0("../../","docs/Source data/","Fig4E",".xlsx"))
pHR_main <- plot_grid(plotlist = pList[names(pList) %in% mainDrug], ncol=1)
pHR_supp <- pList[!names(pList) %in% mainDrug]

For Supplementary Figure

cowplot::plot_grid(plotlist= pHR_supp, ncol=2)

6 Assemble figures for the mansucript

6.1 Main figure 4

pPie <- ggdraw() +  draw_image("../../docs/validation_PieChart_Figure.png", scale = 1) 
pFig4 <- plot_grid(plot_grid(plot_grid(pPie, pVolcano, ncol=1, rel_heights = c(1,0.8), labels = c("A","B"), label_size = 20),
                   pBoxELN_main, ncol=2, rel_widths = c(0.6,1), labels = c("","C"), label_size = 20),
                   NULL,
                   plot_grid(pKM_main, pHR_main, ncol=2, rel_widths = c(1,0.7), labels = c("D","E"), label_size = 20),
                   ncol=1, rel_heights = c(1,0.05,0.6))+
  patchwork::plot_annotation(title = "Figure 4") &  theme(title = element_text(size=25))


writexl::write_xlsx(rownames_to_column(as.data.frame(pVolcano$data),var="rowname"), path=paste0("../../","docs/Source data/","Fig4B",".xlsx"))

pFig4

ggsave("../../docs/Figure_4.pdf", height = 18, width = 14)

6.2 Supplementary plot showing the boxplot and dose-response plot

plot_grid(pBox, pDose, ncol=2, labels=c("A","B"), label_size = 20)

ggsave("../../docs/Figure_ED8.pdf", height = 12, width = 20)

6.3 Supplementary plot showing p-value heatmap and boxplot stratified by ELNrisk

plot_grid(pPheatmap, pBoxELN_supp, ncol=1, labels=c("A","B"), label_size = 20, rel_heights = c(0.5,1))

ggsave("../../docs/Figure_ED9.pdf", height = 20, width =10)

writexl::write_xlsx(rownames_to_column(as.data.frame(pPheatmap$data),var="rowname"), path=paste0("../../","docs/Source data/","ED9A",".xlsx"))

```